Задача этого проекта — научиться предсказывать количество поездок в ближайшие часы в каждом районе Нью-Йорка. Для того, чтобы её решить, сырые данные необходимо агрегировать по часам и районам. Агрегированные данные будут представлять собой почасовые временные ряды с количествами поездок из каждого района. Похожие задачи возникают на практике, если вам необходимо спрогнозировать продажи большого количества товаров в большом количестве магазинов, объём снятия денег в сети банкоматов, посещаемость разных страниц сайта и т.д.

Неделя 2: Работа с геоданными

0. Вступление и подготовка к работе

Авторы проекта предложили использовать инструмент на выбор. Т.к. для меня это второй проект, то я решил его сделать на R. Ввиду того, что в специализации используется Python, я постарался снабдить код подробными комментариями, чтобы (при желании) было легче разобраться.

Для работы потребуются следующие библиотеки:

library(leaflet) # интерактивные карты
library(scales) # в помощь графикам
library(ggplot2) # графики
library(ggmap) # геоданные с Гугла
library(magrittr) # c'est ne pas une pipe
library(data.table) # работа с таблицами данных

# Sys.setlocale('LC_ALL','utf-8') # если неверно отображается кириллица

Информация о версиях библиотек и системе:

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] data.table_1.10.1  magrittr_1.5       ggmap_2.7         
## [4] ggplot2_2.2.0.9000 scales_0.4.1       leaflet_1.0.1     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.8       knitr_1.15.1      maps_3.1.1       
##  [4] munsell_0.4.3     lattice_0.20-34   geosphere_1.5-5  
##  [7] colorspace_1.3-1  R6_2.2.0          rjson_0.2.15     
## [10] jpeg_0.1-8        dplyr_0.5.0       stringr_1.1.0    
## [13] plyr_1.8.4        tools_3.3.2       grid_3.3.2       
## [16] gtable_0.2.0      png_0.1-7         DBI_0.5-1        
## [19] htmltools_0.3.5   yaml_2.1.14       lazyeval_0.2.0   
## [22] rprojroot_1.1     digest_0.6.10     assertthat_0.1   
## [25] tibble_1.2        reshape2_1.4.2    mapproj_1.2-4    
## [28] bitops_1.0-6      htmlwidgets_0.8   evaluate_0.10    
## [31] rmarkdown_1.2     sp_1.2-3          stringi_1.1.2    
## [34] RgoogleMaps_1.4.1 backports_1.0.4   proto_1.0.0

1. Загрузите агрегированные данные о поездках в мае 2016

Просуммируйте общее количество поездок такси из каждой географической зоны и посчитайте количество ячеек, из которых в мае не было совершено ни одной поездки.

# загрузка
regions <- fread('data_in/regions.csv') # данные о регионах
dt_agg <- fread('data_out/2016_05.csv')

# суммирование по регионам
dt_month_regions <- dt_agg[, .(n = sum(n)), by = .(region)]

dt_month_regions[n == 0, .N]
## [1] 1283

Таким образом из 1283 регионов не было совершено ни одной поездки.

##########################
# Константа с геоданными #
##########################

# Нью-Йорк вписан в прямоугольник от -74.25559 до -73.70001 градусов долготы и от 40.49612 до 40.91553 широты
NY <- list(lon = c(-74.25559, -73.70001), lat = c(40.49612, 40.91553))
NY$CELLS <- 50 # число ячеек по широте или долготе
NY$lon_step <- (NY$lon[2] - NY$lon[1]) / NY$CELLS # шаг по долготе
NY$lat_step <- (NY$lat[2] - NY$lat[1]) / NY$CELLS # шаг по широте
NY$lon_breaks <- seq(NY$lon[1], NY$lon[2], by = NY$lon_step) # границы регионов по долготе
NY$lat_breaks <- seq(NY$lat[1], NY$lat[2], by = NY$lat_step) # границы регионов по широте

2. Нарисуйте статическую карту Нью-Йорка

Поставьте на карте точку там, где находится Эмпайр-Стейт-Билдинг.

Для отображения статических карт будем пользоваться пакетами ggmap и ggplot2.

# координаты Empire State Building
# geocode('Empire State Building')

# введем координаты вручную
ESB <- list(lon = -73.98566, lat = 40.74844)

# скачиваем карту с маркером Эмпайр-Стейт-Билдинг
NY$map <- get_googlemap(c((NY$lon[1] + NY$lon[2]) / 2, 
                          (NY$lat[1] + NY$lat[2]) / 2),
                        zoom = 10, scale = 2,
                        size = c(640, 640),
                        markers = data.frame(lon = ESB$lon, lat = ESB$lat)
                        )

# отображаем и обрезаем карту по границам регионов
ggmap(NY$map, extent = 'device') +
    scale_x_continuous(limits = c(NY$lon[1], NY$lon[2])) +
    scale_y_continuous(limits = c(NY$lat[1], NY$lat[2])) + 
    labs(title = 'Карта Нью-Йорка', 
         subtitle = 'Эмпайр-Стейт-Билдинг отмечен маркером')

3. Визуализируйте данные о поездках

Поверх статической карты Нью-Йорка визуализируйте данные о поездках из каждой ячейки так, чтобы цветовая шкала, в которую вы окрашиваете каждую ячейку, показывала суммарное количество поездок такси из неё.

Для начала посмотрим на распределение общего числа поездок по регионам за месяц. Для наглядности отобразим квадратный корень и логарифм этой величины:

# гистограмма n
dt_month_regions[n>0, .(region, n, sqrt_n = sqrt(n), log_n = log(n+1))] %>%
    melt(measure.vars = 2:4, id.vars = 1) %>%
    ggplot(aes(value, ..ncount..)) + 
    geom_histogram(bins = 15) + 
    facet_wrap(~variable, scales = 'free_x') + 
    labs(title = 'Гистограммы числа поездок', 
         subtitle = 'n, sqrt(n) и log(n+1)',
         x = 'Значение', y = 'Норм. частота') + 
    scale_y_continuous(breaks = c(0, 0.5, 1))

Очевидно, что если не преобразовывать число поездок, то практически вся карта будет раскрашена одноцветно, т.к. много нулей и маленьких значений. Поэтому для раскраски каждого региона воспользуемся логарифмическим преобразованием, а чтобы ячейки с небольшими значениями имели меньший цветовой вес, пусть прозрачность определяется корнем величины. Регионы с нулевым значением отрисовывать не будем. В этом случае, тепловая карта суммарного числав поездок будет выглядеть так:

# добавляем каждому региону его границы
dt_month_regions <- cbind(dt_month_regions, regions[, 2:5])

# скачиваем чистую карту
NY$map <- get_googlemap(c((NY$lon[1] + NY$lon[2]) / 2, 
                          (NY$lat[1] + NY$lat[2]) / 2),
                        zoom = 10, scale = 2,
                        size = c(640, 640)
                        )

# отображаем карту
ggmap(NY$map, extent = 'device') + 
    # рисуем квадратики, исключаем нулевые
    geom_rect(data = dt_month_regions[n > 0], 
              aes(xmin = west, xmax = east, ymin = south, ymax = north,
                  fill = log(n), alpha = sqrt(n)), 
              size = 0.1, color = 'gray60', inherit.aes = FALSE) +
    # обрезаем границы
    scale_x_continuous(limits = c(NY$lon[1], NY$lon[2])) +
    scale_y_continuous(limits = c(NY$lat[1], NY$lat[2])) +
    # заливка квадратиков
    scale_fill_distiller(palette = "OrRd", breaks = pretty_breaks(n = 6),
                         direction = 1) +
    # прозрачность квадратиков в зависимости от n
    scale_alpha(range=c(0.2, 0.5)) +
    # убираем легенды
    theme(legend.position="none") +
    # заголовок
    labs(title = 'Тепловая карта поездок за месяц')

На карту попали регионы, откуда поездок быть не может (бухты, река).

4. Вставьте интерактивную карту Нью-Йорка

Вставьте интерактивную карту Нью-Йорка — такую, которую можно прокручивать и увеличивать. Поставьте метку там, где находится статуя свободы.

Для интерактивных карт воспользуемся пакетом leaflet.

# координаты Empire State Building
# geocode('Statue of Liberty')

# введем координаты вручную
SoL <- data.frame(lon = -74.0445, lat = 40.68925, name = 'Statue of Liberty')

# отображаем карту
leaflet(SoL, width = '100%') %>% addTiles() %>% 
    addMarkers(~lon, ~lat, popup = ~name) %>%
    setView(SoL$lon, SoL$lat, zoom = 12)

5. Отобразите среднее за месяц количество поездок

Нарисуйте на интерактивной карте Нью-Йорка ячейки так, чтобы их цвет показывал среднее за месяц количество поездок такси в час из этой зоны.

Прежде всего необходимо рассчитать среднее количество поездок в час для каждого региона. Попробуем отобразить “горячую карту” так же, как и статичную карту.

Рассчитаем прозрачность как корень из числа поездок и отмасштабируем множество этих значений на отрезок [0.2; 0.6].

Для получения цвета для каждой зоны, создадим палитру из 50 оттенков (как в предыдущей тепловой карте). Среднее число поездок из каждого региона прологарифмируем, а все множество логарифмов разобьем на 50 интервалов, сопоставив таким образом каждому логарифму среднего числа поездок свой оттенок из палитры. Для удобства добавим текст по клику для квадратиков, а для самой карты предусмотрим возможность отключения как слоя. По аналогии с предыдущей картой, не будем отображать регионы, среднее значение которых меньше 0.1.

# рассчитываем среднее количество поездок
dt_agg_avr <- dt_agg[, .(avr_n = mean(n)), by = .(region)]
# добавляем границы квадратиков
dt_agg_avr <- cbind(dt_agg_avr, regions[,2:5])


# палитра для квадратиков 
shades <- 50 # 50 оттенков красно-оранжевого
pal_OrRd <- seq(0, 1, length.out = 50) %>% 
    (brewer_pal(palette = 'OrRd', direction = -1)(9)[1:5] %>% 
         rev %>% gradient_n_pal)

# прологарифмируем средние значения и разобьем их на 50 градаций 
# сопоставляем цвет из палитры каждому полученному числу
dt_agg_avr <- dt_agg_avr[, 
    fill_color := pal_OrRd[log(avr_n + 1) %>% cut(breaks = shades) %>%
                               as.numeric()]]

# добавляем признак прозрачности, ~ sqrt(n)
dt_agg_avr <- dt_agg_avr[, alpha := sqrt(avr_n + 1) %>% rescale(c(0.2, 0.6))]

# отрисуем карту
leaflet(dt_agg_avr[avr_n > 0.1], width = '100%') %>% addTiles() %>% 
    # квадратики
    addRectangles(lng1 = ~west, lng2 = ~east, lat = ~south, lat2 = ~north,
                  # границы квадратиков
                  weight = 1, color = 'gray', opacity = 4/10,
                  # заливка квадратиков
                  fillOpacity = ~alpha, fillColor = ~fill_color,
                  # описание квадратика
                  popup = ~paste('region:', region, 'avr_n:', round(avr_n, 1)),
                  # для управления слоями
                  group = 'heatmap') %>%
    # чтобы квадратики можно было убирать
    addLayersControl(
        overlayGroups = c("heatmap"),
        options = layersControlOptions(collapsed = FALSE))

6. Отфильтруйте ячейки

Чтобы не выбирать из всех 2500 ячеек вручную, отфильтруйте ячейки, из которых в мае совершается в среднем меньше 5 поездок в час. Посчитайте количество оставшихся. Проверьте на карте, что среди этих ячеек нет таких, из которых поездки на самом деле невозможны.

После того, как уберем зоны, откуда было в среднем было меньше 5 поездок в час, останется 102 региона:

dt_agg_avr[avr_n >= 5, .N]
## [1] 102

Отобразим и выделим их на карте. Включая и отключая слой с регионами можно убедиться, что из всех этих зон возможны поездки.

# отрисуем карту
leaflet(dt_agg_avr[avr_n >= 5], width = '100%') %>% addTiles() %>% 
    # квадратики
    addRectangles(lng1 = ~west, lng2 = ~east, lat = ~south, lat2 = ~north,
                  # границы квадратиков
                  weight = 1, color = 'gray', opacity = 3/4,
                  # заливка квадратиков
                  fillOpacity = 1/2, fillColor = 'maroon',
                  # описание квадратика
                  popup = ~paste('region:', region, 'avr_n:', round(avr_n, 1)),
                  # для управления слоями
                  group = 'heatmap') %>%
    # чтобы квадратики можно было убирать
    addLayersControl(
        overlayGroups = c("heatmap"),
        options = layersControlOptions(collapsed = FALSE))

7. Опубликуйте ноутбук

Ноутбук опубликован на rstudioconnect.com. Файлы проекта также можно найти и на гитхабе.